#include "reconstruction.hpp"

namespace hms
{

/**
 * @brief Combines Murillo's condition (d_min >= d_max - d_min)
 * and a check for dry cells on either side. Returns true if both are passed,
 * meaning that central or high order reconstruction is possible,
 * and false otherwise, meaning that first order reconstruction is required.
 * 
 * @tparam SolverType required for SolverType::minWaterDepth
 * @param dCellLR 
 * @return true 
 * @return false 
 */
template<typename SolverType>
bool murillo( const Array2s& dCellLR ){
	scalar
		dMin { dCellLR.minCoeff() },
		dMax { dCellLR.maxCoeff() };

	return (
		dMin >= SolverType::minWaterDepth &&
		dMin >= dMax - dMin
	);
}

template<EdgeRange orien, typename MeshType>
std::array<scalar, 2> dxQuotient(
	const MeshBase<MeshType>& mesh,
	[[maybe_unused]] Index cellLL
){
	using ER = EdgeRange;

	if constexpr ( std::is_same_v<MeshType, UniMesh> ){
		return {{ 1, 1 }};
	} else if constexpr ( std::is_same_v<MeshType, RectMesh> ){
		RectCell llCell { mesh.cell( cellLL ) };
		constexpr Index ng {1};
		Array3s dx;
		if constexpr ( orien == ER::horizontal ){
			/* NEEDS TESTING */
			dx = mesh.cells().dy().segment( llCell.j() + ng, 3 );
		} else {
			/* NEEDS TESTING */
			dx = mesh.cells().dx().segment( llCell.i() + ng, 3 );
		}
		return {{
			dx[1] / dx[0],
			dx[1] / dx[2]
		}};
	} else {
		assert(false);
		return {{ 0, 0 }};
	}
}

template<typename MeshType>
template<ReconstructionType recType, EdgeRange orien>
EdgeValues SWESolver<MeshType>::reconstruction(
	const MeshBase<MeshType>& mesh,
	const typename MeshType::EdgeType& edge,
	const Array2i& lr,
	const RowArray2s& zCellLR,
	const RowArray2s& dCellLR
){
	using RT = ReconstructionType;
	using S = SWESolver<MeshType>;
	static constexpr bool
		isUniMesh  { std::is_same_v<MeshType, UniMesh> },
		isRectMesh { std::is_same_v<MeshType, RectMesh> },
		isOrtho    { isUniMesh || isRectMesh };

	auto getFieldVals = [&](const auto& field, const Array2i& cells){
		return field.values()(hms::all, cells);
	};

	ArrayNNs<2,2> vCellLR { getFieldVals( v(), lr ) };

	EdgeValues ev;

	using namespace reconstruction;

	if ( murillo<S>(dCellLR) && isOrtho && recType != RT::firstOrder ){

		RowArray2s hCellLR { zCellLR + dCellLR };
		ArrayNNs<2,2> qCellLR { vCellLR.rowwise() * dCellLR };

		if constexpr (recType == RT::secondOrder){
			Array2i llrr { edge.cells(2) };
			RowArray2s
				dCellLLRR { getFieldVals(d(), llrr) },
				zCellLLRR { getFieldVals(z(), llrr) },
				hCellLLRR { zCellLLRR + dCellLLRR };

			ArrayNNs<2,2>
				qCellLLRR { getFieldVals(v(), llrr).rowwise() * dCellLLRR };

			/* scalar */
			auto [dxQuotL, dxQuotR] = dxQuotient<orien>(mesh, llrr[0]);

			reconstruction::highOrder<S>(
				ev,
				dCellLR, dCellLLRR, hCellLR, hCellLLRR, qCellLR, qCellLLRR,
				dxQuotL, dxQuotR, this->settings().limiterFunction
			);
		} else if constexpr (recType == RT::central) {
			reconstruction::central<S>(ev, dCellLR, hCellLR, qCellLR);
		}
	} else {
		reconstruction::firstOrder<S>( ev, dCellLR, zCellLR, vCellLR );
	}
	reconstruction::hydrostatic_firstOrder(ev);

	return ev;
}

} // namespace hms